Random-exchange quantum Heisenberg antiferromagnets on the square lattice 
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The ground state properties of random-exchange spin- 1/2 Heisenberg antiferromagnets on the 
square lattice are investigated using a combination of quantum Monte Carlo simulations, exact 
numerical diagonalizations, and spin wave calculations. Whereas arbirarily weak disorder has a 
drastic effect on Id Heisenberg AFM, we find that in two dimensions the characteristics of the ground 
state like long-range order is robust even against strong disorder. While the antiferromagnetic order 
parameter and the spin stiffness are exponentially reduced for singular exchange distributions, they 
vanish only in the limit of infinite randomness. 

PACS numbers: 75.40.Mg, 75.10.Nr, 02.70.Ss, 75.30.Ds 



The spin-1/2 Heisenberg antiferromagnet (AFM) on 
the square lattice has attracted a lot of interest in the 
last two decades [H S H S 0, 0, 0, motivated in 
particular by the suggestion that the high-temperature 
superconductivity in cuprates is related to the magnetic 
properties of their CuO planes 0]. Its ground state (GS) 
is antiferromagnetically ordered and the value for the 
staggered magnetization jtiaf = 0.3070(3) is known to 
high accuracy Recently it was found that this mag- 
netic long range order (LRO) is robust against the intro- 
duction of static non-magnetic impurities, as observed in 
Mg or Zn doped La2Cu04 0- In fact, using numeri- 
cal simulations, it was shown that LRO persists up to 
the percolation threshold for site dilution [HI . Fur- 
thermore, the critical exponents for the transition to the 
paramagnetic phase are those of the classical percolation 
transition [l^ . 

The question then arises, how the GS properties 
change for other forms of quenched disorder, in partic- 
ular for generic bond disorder, i.e. randomness in the 
strength of the antiferromagnetic exchange coupling. In 
the one-dimensional (ID) case, where the GS of the pure 
system is critical, bond disorder is relevant [l^ . and an 
infinitesimal amount of bond disorder drives the system 
into the random sing let phase with unconventional scal- 
ing properties 0,0,0] . On the other hand, the gapped 
spin-liquid GS of antiferromagnetic spin-1/2 Heisenberg 
ladders is stable against the introduction of bond disor- 
der 0, 01 . Indications for similar stability were ob- 
tained recently also for two-dimensional (2D) Heisenberg 
antiferromagnets using the strong disorder rcnormaliza- 
tion group |19|,l2^], which is, however, reliable only for 
the opposite scenario, i.e. when the GS is described by 
an infinite randomness fixed point scenario. 

Here, we employ exact numerical methods to reliably 
estimate the stability of the GS of 2D Heisenberg AFMs 
against exchange-randomness. In particular, we consider 
the bond-disordered spin-1/2 Heisenberg AFM on the 



square lattice, defined by the Hamiltonian 

H = ^ Jij Si ■ Sj, (1) 

where {i, j) are nearest neighbor bonds on a L x L square 
lattice {N — L^), Si are spin-1/2 operators, and the 
AFM exchange couplings > are quenched random 
variables with probability distribution 'P(J). In the fol- 
lowing, we will consider two different types of disorder 
distributions: (i) a VF-dependent flat, bounded bond dis- 
tribution, where the bonds are uniformly and symmetri- 
cally distributed around 1 with width W: 

v{j) = e( J - (1 - w)) ■ e((i + w)- J) 1 2w (2) 

and (ii) a singular distribution controlled by the param- 
eter 5 and given by 

v{j) = j-^+^''5-^ ■ e(j)e(i - J) . (3) 

This second type of disorder distribution allows us to 
study the stability of the GS towards a singular prolifer- 
ation of weak bonds in the lattice structure. 

We analyze the GS properties of the model using 
a combination of quantum Monte Carlo (QMC) simula- 
tions, exact diagonalizations (ED) based on the Lanczos 
method, and spin wave theory (SW). In particular, we 
consider the stability of the staggered magnetization and 
the spin stiffness towards bond disorder, using disorder 
averaged quantities as obtained from several thousand 
disorder realizations from the distributions (2) and (3) 
with different parameters W or 5, respectively. 

A Nccl-ordered GS in the thermodynamic limit is sig- 
naled by a non-vanishing infinite system size limit TOaf 
of the staggered structure factor per site s(7r, tt): 

= ^ {<^yiySrf'^ ^ mip {N oo). 

(4) 
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FIG. 1: (Color online) Disorder averaged value of the GS 
staggered structure factor (s(7r, tt)) as a function of 
Filled symbols are ED data and open symbols QMC data. 
The number of random samples varied between 1000 for the 
largest lattice and 10000 for the smaller lattices. Dashed lines 
denote third order polynomial fits of the finite size data. The 
inset exhibits the convergence of (s(7r,7r)) using /3-doubling 
scheme for 5 = 1, 2, 3, 4, 5 (top to bottom) and L = 32. 



The QMC simulations have been performed using the 
stochastic series expansion technique vifith a directed loop 
update scheme jl^. For the clean (i.e. no disorder) case, 
we studied square lattices up to L = 64. We found that 
a scaling of the inverse temperature as /? = 8 x L allows 
to obtain GS properties. For the random case we needed 
to perform the QMC simulations at much lower temper- 
ature, as correlations on weak bonds (i.e. small Jij) still 
develop at very low temperatures. To stimate GS quan- 
tities reliable, we used the /3-doubling scheme 0j which 
accelerates the simulations considerably at low temper- 
atures. We were able to extract the GS properties of 
systems with up to 32 x 32 lattice sites, for values of the 
disorder-parameter up to 5 = 5 by doubling /3 up to 2^^ 
(see the inset of Fig. For the small lattices we have 
also carefully checked the convergence against the ED 
results. The number of random samples considered for 
disorder averaged quantities varied between 1000 for the 
largest lattice and 10000 for smaller lattices. 

In Fig. n we show the disorder averaged GS values of 
(s(7r, tt)), which we obtained from the QMC and ED cal- 
culations as functions of 1/L. For the clean case (8 = 0) 
QMC data for lattices ranging from L = 4 up to L = 64 
have been extrapolated to the thermodynamic limit us- 
ing an unconstrained third order polynomial fit in 1/L. 
The obtained staggered magnetization ttzaf = 0.3064(2), 
is in good agreement with previous estimates |^ [i^ . In- 
creasing the randomness (5, we found that large system 
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FIG. 2: (Color online) Results from spin wave theory for the 
disorder averaged value of the GS staggered structure factor 
(s(7r, tt)) as a function of 1/L. Dashed lines denote polynomial 
fits in 1/L. 



sizes are required to reliably extrapolate to the thermo- 
dynamic limit. While in the clean case and for weak 
disorder the finite size scaling is dominated by a term 
linear in 1/L, higher-order terms become important for 
larger values of 8. This can be seen from Fig.^ where for 
strong disorder and large systems sizes, the curvature of 
the polynomial extrapolations become more pronounced, 
indeed requiring up to third order polynomials in 1/L. 
We checked that the extrapolated values did not change 
within the statistical error bars upon allowing higher or- 
der fitting polynomials. 

To further analyze the effects of the randomness on 
the quantum fiuctuations above the classical Neel state, 
we employ linear spin-wave |3| theory. This appears 
justified given the finite, although strongly reduced, stag- 
gered moment obtained within the QMC simulations. As 
translational symmetry is broken in each realization for 
the disordered system, we adopt a a real-space formu- 
lation of the spin- wave approach to finite systems. In 
the following, we use the method of Ref. |25j, which is 
based upon solving a non-hermitian eigenvalue problem 
and subsequent orthogonalization to obtain the bosonic 
eigenmodes of the quadratic part of the spin- wave Hamil- 
tonian. Similar approaches have been used recently to 
study site-dilution effects in the Heisenberg model on the 
square [2^ and honeycomb [i^l lattice. For details on the 
numerical scheme, we refer to Ref. |25| . 

We perform the above procedure for typically 1000 dis- 
order realizations for system sizes L = 8, 16, 24, 32 and 
40 for various values of 8. The spin- wave estimate of 
(s(7r,7r)) for given values of L and b is finally obtained 
by averaging the local staggered moments from the vari- 
ous disorder realizations. We found that for (5 > 1 a small 
fraction of sites attain negative values of the local stag- 
gered moment. We treat this artifact of the linear spin- 
wave theory, related to the singular disorder distribution 



0.4 



~<0.2 



0.1 



• QMC 

-- ~exp(-0.3 5) 

o spin-waves 

-- ~exp(-0.26 8) 




1 



FIG. 3: (Color online) Disorder averaged antiferromag- 
netic order parameter as a function of the disorder strength 
S. Dashed lines denote exponential fits of the data for 5 > 1, 
and the inset shows the same data on a semi- logarithmic plot. 
Note also that error bars are smaller than symbol sizes. 



at J = 0, by explicitly setting the local staggered mo- 
ments of the corresponding sites to zero before perform- 
ing the averaging procedure. Results of the SW calcula- 
tion are shown in Fig.|21for various disorder strengths as 
functions of l/L. Extrapolations to the thermodynamic 
limit have been performed using third order polynomial 
fits in l/L, where again high-order terms become impor- 
tant as the randomness increases. 

Based on our finite size scaling analysis, we now con- 
sider the behavior of the staggered magnetization tuaf 
in the presence of bond-randomness. We plot the val- 
ues obtained from both the ED and QMC and the SW 
calculations as functions of S in Fig. O 

We find the Neel order is very robust against bounded 
bond randomness. While the SW approach slightly over- 
estimates the order parameter (as expected), the form 
of the 5-dependence of mAp calculated within SW the- 
ory agrees well with the ED and QMC results. Using 
the distribution of Eq. (3), the order parameter shows an 
exponential decay (c.f. the inset Fig. O of the form 



niAF ~ exp(-c,„J), 



(5) 



with c,n = 0.301(5) obtained from the ED and QMC 
calculations, and Cm = 0.260(5) within SW theory. This 
suggests that the robust LRO vanishes only in the limit 
of infinite randomness, i.e. for 5 — + oo. 

In addition to the staggered magnetization, the or- 
dered nature of the Heisenberg antiferromagnet is accom- 
panied by a finite value of the spin stiffness, defined as 
the second derivative of the GS energy with respect to a 
twist angle $ introduced at the boundary along a direc- 
tion perpendicular to the order parameter: 
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FIG. 4: (Color online) Disorder averaged values of the spin 
stiffness p7 obtained from QMC simulations on systems with 
up to L = 32 and carefuU finite size scaling to the thermody- 
namic limit (see text). The solid line denotes an exponential 
fit Eq. (7). Inset; Finite size scaling of the disorder average 
LJ — current correlator A2{L) = — eo(Z/)/3 — 2ps{L)/3 plot- 
ted vs l/L for 5 = 0, 1, 2, 3, 4, 5 (different symbols from 
top to bottom). Dashed lines are quadratic fits. Error bars 
are smaller than symbol sizes. 



In contrast to mAp, the spin stiffness is a dimensionful 
quantity, and in the clean case, where ps = 0.175(2)J, 
scales proportional to the exchange constant J. A naive 
expectation would suggest that in the disordered case, ps 
scales with the averaged value of the exchange, J ~ 
leading to an algebraic decay for large values of S. In the 
following, we show that instead ps shows an exponen- 
tial decay, similar to the staggered magnetization, toap. 
From the QMC simulations, ps can be evaluated for each 
disorder realization by measuring the world-line winding 
number fluctuation |3| , and finally performing a disorder 
averaging. In a clean system, a twist $ enforced at the 
boundaries will be equally distributed among the bonds 
along the corresponding direction, resulting in a homo- 
geneous twist in the order parameter of <I>/L per bond. 
In contrast, for a bond-disordered system, larger twists 
occur for weaker bonds |28l |. resulting in an inhomoge- 
neous distribution of the local twist angle. The winding 
number however, being a global quantity, provides an es- 
timate for the global rigidity of the system under such a 
twist in the boundary conditions. 

The spin stiffness can be expressed as = — ^ ~ , 
in terms of GS energy per site without twist, eo, and A^, 
the Lo = current-current correlation function As(cj) = 

Jo dTe-^^^jsiTYLm H 113. According to spin- 
wave calculations |22|,|21|, in the clean case A°(i) scales 
like AO(L) = AO -h ^ fl + . . .^hereas the GS en- 
ergy scales as cq (L) — cq ^ ^ H, 0| . We find that the 
scaling of eo{L) also holds for the disordered aver- 
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aged values. To reliably obtain the disorder average spin 
stiffness pj in the thermodynamic limit, we subtracted 
the energy contribution to the stiffness and performed 
polynomial fits in 1 / L for the disorder averaged values of 
AO(L), as shown in the inset of Fig. 01 Interestingly, for 
increasing disorder, finite size effects get strongly reduced 
which implies that the thermodynamic limit value for p5 
is, within error bars, reached for smaller systems as disor- 
der increases. Since more QMC steps are needed for each 
random sample to compute the fluctuations in the wind- 
ing numbers, than for the structure factor, the disorder 
average was performed over slightly less random samples. 
Typically 1000 samples were used for the smaller and 200 
for the largest lattice with L = 32. Our results for pj as 
a function of 5 arc shown in Fig. 01 Wc find that the 
spin stiffness remains finite upon increasing the disorder, 
showing an exponential reduction with J, similar to the 
AF order parameter. In fact, a decay 

/9s(<5) exp(-Cp(5), (7) 

with Cp = 1.60(2), fits the QMC data for 5 > I. 

In conclusion we studied the ground state properties 
of the bond-disordered spin- 1/2 Heiscnberg model on 
the square lattice, and found the antifcrromagnctic or- 
der to be rather robust against bond-randomness: Both 
the staggered magnetization and the spin stiffness get re- 
duced exponentially with 5, but vanish only in the limit 
of infinite randomness. These observations suggest a re- 
duced relevance of disorder effects on quantum Heisen- 
berg models in dimensions D > 1. Indeed, and in con- 
trast to the Ising model in a transverse field, where dis- 
order is relevant also for D > 1 j^^j the infinite 
randomness fixed point is unstable in XY and Heiscn- 
berg AFMs for D > 2 0. This result of a rcnormal- 
ization group study is confirmed here for the Heiscnberg 
case on the square lattice. It would be interesting to 
achieve an analytical understanding of the exponential 
decay wc obtained for both the order parameter and 
the spin stiffness. In an earlier study, performed for a 
similar random-exchange AFM, a finite critical disorder 
strength was found to destroy Neel order and to drive 
the system into a paramagnetic phase [l^l- Based on 
simulations using significantly larger system sizes at es- 
sentially zero temperature, we conclude on the absence 
of any such finite critical disorder strength. We further- 
more expect a similar robustness of Neel-ordered GSs of 
quantum AFMs on other bipartite lattices, such as the 
honeycomb lattice. In contrast, for the case of a trian- 
gular lattice the interplay between randomness and frus- 
tration is expected to destroy the ordered GS of the pure 
case [s^ , and drive the system into a spin-glass phase for 
strong bond-disorder [l9j. Finally, the question arises, 
how the excitation spectrum and finite temperature prop- 
erties are affected by bond-disorder. Preliminary results 
|36l | indicate that upon increasing the disorder strength, 
the character of the low-energy excitations changes from 



extended spin waves to localized spin-cluster flips. 
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